function fsw_proto()
  clear 
  global N k x diagnostics y0 g L y0_n
  format long e
  cfl = 0.25;
  N = 512;  
  L = 20*pi;
  c = 2*pi/L;
  dt = cfl*(L/N);
  tmax = 200; 
  nsteps = floor(tmax/dt)+1;
  dt = tmax/nsteps;
  %%%%%%%%%%%%%%%%%%%%%%%            Parameters
  A = 0.1;   % elevation 
  B = 0.000;   % velocity
  k0 = 2;
  y0 = -5; 
  y0_n = y0;
  y0_s = y0;
  g = 1;
  kp = 6;
  kc = 100*c;
  diagnostics = 0;
  evolve = 1;
  k = c*(-N/2:N/2-1);
  x = L*(k/c)/N;
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  z_k = A*exp(-(k/k0).^2).*exp(0.*2i*pi*rand(1,N)).*exp(-1i*k.*L/3);
  %z_k = zeros(N,1); z_k(N/2+2) = A;
  a_k = B*exp(-(k/k0).^2).*exp(0.*2i*pi*rand(1,N)).*exp(1i*k.*L/3);
  
  y_k = herm(z_k)';
  p_k = herm(a_k)';
  z_k = 0+1i*y_k;
  
  y_k(N/2+1) = y0;
  p_k(N/2+1) = 0.*real(p_k(N/2)+p_k(N/2+2));
  z_k(N/2+1) = 1i*y0;
  
  
  p_k = 1.*p_k;
  
  v_y0 = zeros(N,1);
  z = b_ft(z_k,'nonsymmetric');
  y = b_ft(y_k,'symmetric');
  p = b_ft(p_k,'symmetric');
  fid = fopen('harmonics_10.txt','w');  
  if evolve == 1
   t = 0;
   q(:,1) = y_k(:);
   q(:,2) = p_k(:);
   q(:,3) = z_k(:);
   q(:,4) = p_k(:);
   for j = 1:nsteps
     rq = rhs_linear(q);
     qt = q + 0.5*dt*rq;
     rq = rhs_linear(qt);
     q = q + dt*rq;     
     
     q(:,3) = q(:,3).*exp(-dt*(k(:)/kc).^kp);
     q(:,4) = q(:,4).*exp(-dt*(k(:)/kc).^kp);

     %y0  = q(N/2+1,1);
     y0  = imag(q(N/2+1,3));
     t = t + dt;

     v_y0(N/2+1) = y0; 
     
     
     %zt = b_ft(q(:,3),'nonsymmetric');
     xt_n = b_ft(herm(q(:,3)),'symmetric');
     dxt_n = b_ft(herm(1i.*k(:).*q(:,3)),'symmetric');

     ml = trapz((b_ft(antiherm(q(:,3)),'symmetric') ).*(1+xt_n(:)))*(L/N);
     ml_k = sum(antiherm(q(:,3)).*( v_y0/y0 + 1i.*k(:).*herm(q(:,3)) )  )*L;

     Kin = 0.5*sum(k(:).*tanh(-y0*k(:)).*abs(q(:,4)).^2)*L;
     %dxt = b_ft(1i*k(:).*That(q(:,1)),'symmetric');
     
     
     
     Pot = trapz(0.5*g*(b_ft(antiherm(q(:,3)),'symmetric')).^2.*(1+dxt_n))*(L/N);
     %Pot_2 = sum(0.5*g*(q(:,1)-v_y0).*conj(q(:,1)-v_y0))*L;
     %ml = sum(conj(antiherm(q(:,3))).*( 1i*k(:).*That(q(:,1)) )   )*L;
     %fprintf('T = %e\tKin = %e\tPot = %e\tTotal = %e\t Total2 = %e\n', t, Kin, Pot, Kin + Pot, Kin + Pot_2)
     fprintf('T = %e\tdel y0 = %e\tdepth = %e\tK = %e\tP = %e\tT = %e\tml = %e\n', ...
              t, y0-y0_s, -y0, Kin, Pot, Kin+Pot, ml)
     fprintf(fid,'%e\t%e\t%e\t%e\t%e\t%e\n', t, Kin, Pot, Kin + Pot,ml, ml_k);
     
     figure(1)
     semilogy(k(:), [ abs(q(:,3)) abs(q(:,4)) ]);
     axis([k(1) k(end) 1e-22 1])
     
     figure(2)
     plot(x(:),  real(b_ft(q(:,1),'symmetric'))-y0,...
          x(:)+xt_n(:),  imag(b_ft(q(:,3),'nonsymmetric'))-y0, ...
          x(:), 0.*b_ft(1i*k(:).*q(:,2),'symmetric')  );
     axis([-L/2 L/2 -10 10])
     

     pause(0.0005) 
   end
  end
  fclose(fid);
  if diagnostics == 1
      
%     rp_k = Rhat(p_k);
%     tp_k = That(p_k);
%   
%     dp = b_ft(dp_k,);
%     rp = b_ft(rp_k);  
%     tp = b_ft(tp_k);
%     pp_k = That(rp_k);
%     pp = b_ft(pp_k);
      
    Kin = 0.5*sum(k(:).*tanh(-y0*k(:)).*abs(p_k(:)).^2)*L;
    
    v_y0(N/2+1) = y0; 
    y = b_ft(y_k(:,1),'symmetric');
    dxt = b_ft(1i*k(:).*That(y_k(:,1)),'symmetric');
    Pot = trapz(0.5*g*(y-y0).^2.*(1+0.*dxt))*(L/N);
    Pot_2 = sum(0.5*g*(y_k-v_y0).*conj(y_k))*L;

    fprintf('Kin = %e\tPot = %e\tPot k = %e\tTotal = %e\t Total2 = %e\n', Kin, Pot, Pot_2, Kin + Pot, Kin + Pot_2)

    
    
    
    figure(1)
    %plot (x, [real(z); imag(z); real(p)'; imag(p)'] )
    plot (x+0.*dxt', [real(y)'; imag(y)'; real(p)'; imag(p)'] )

    axis([-L/2 L/2 -2 2])
%    axis([-L/2 L/2 y0 0.5])
    legend('rz','iz','rp','ip')
    figure(2)
    semilogy (k, ([abs(real(y_k))'; abs(imag(y_k))'; abs(real(p_k))'; abs(imag(p_k))']) )
    legend('ry_k','iy_k','rp_k','ip_k')
    %plot (k, [real(z_k); imag(z_k); real(p_k)'; imag(p_k)'] )
    %legend('rz_k','iz_k','rp_k','ip_k')
    %figure(3)
    %plot(x,[p; pp; 0.*dp; 0.*gradient(p,L/N) ;0.*rp; 0.*tp])
  end  
    
  
  
end
function out = herm(in)
  global N;
  out = zeros(N,1);
  for j = 1:N/2-1
      out(N/2+1+j) = 0.5*(in(N/2+1+j)+conj(in(N/2+1-j)));
      out(N/2+1-j) = conj(out(N/2+1+j));
  end
  out(1) = 0;
end
function out = antiherm(in)
  global N;
  out = zeros(N,1);
  for j = 1:N/2-1
      out(N/2+1+j) = -0.5i*(in(N/2+1+j)-conj(in(N/2+1-j)));
      out(N/2+1-j) = conj(out(N/2+1+j));
  end
  out(1) = 0;
  %out(N/2+1) = 0; except zero harmonic must stay 1i*y0

end
function out = rhs_linear(in)
  global k y0 g N 
  
  dp  = b_ft(1i*k(:).*in(:,4),'symmetric');
  rdp = b_ft(-k(:).*tanh(-y0*k(:)).*in(:,4),'symmetric');
  J = abs(1 + b_ft(1i.*k(:).*in(:,3),'nonsymmetric')).^2;
  ks = f_ft(rdp./J);
%  z_u = b_ft(1i.*k(:).*in(:,3));
  %avg_st = trapz(imag(z_u)./real(z_u))*(L/N);
  
  A = 0.5*(dp.^2 - rdp.^2)./J;
  B = dp.*b_ft(That(ks),'symmetric');
  C = (1+b_ft(1i.*k(:).*in(:,3),'nonsymmetric')).*b_ft(1.i*ks + That(ks),'nonsymmetric');
  
  out(:,1) =  k(:).*tanh(-y0*k(:)).*in(:,2);
  out(:,2) =  -g*in(:,1);
  out(N/2+1,2) = 0; % otherwise zero mode is growing as -y0*t in time ( bad for numerics)
  
  %out(:,3) = -1.i*k(:).*(1-tanh(-y0*k(:))).*in(:,4);
  out(:,3) = -0.*f_ft(C);                              % here
  out(:,4) = -0.*g*antiherm(in(:,3))-f_ft(A(:)+B(:));  % here
  %out(N/2+1,4) = 0; % otherwise zero mode is growing as -y0*t in time ( bad for numerics)
  
end
function out = f_ft(in)
  global N
  out = fftshift(fft(in))/N;
end
function out = b_ft(in,method)
  global N
  out = ifft(ifftshift(in),method)*N;
end
function out = That(in)
  global k y0 N;
  out = zeros(N,1);
  out(:) = -1i.*in(:)./tanh(-y0*k(:));
  out(N/2+1) = 0;
end
function out = Rhat(in)  
  global k y0 N;
  out = zeros(N,1);
  out(:) = 1i.*in(:).*tanh(-y0*k(:));
end
function out = der(in)
  global k N;
  out = zeros(N,1);
  out(:) = 1i.*k(:).*in(:);
end